{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Covariance Propagation\n",
"\n",
"This tutorial demonstrates how to propagate state covariance alongside position and velocity using satkit's `satstate` class. Covariance propagation is essential for:\n",
"\n",
"- **Uncertainty quantification**: Understanding how errors in an initial state estimate grow over time\n",
"- **Conjunction assessment**: Determining the probability of collision between two objects requires knowing the uncertainty in both objects' positions\n",
"- **Sensor tasking**: Deciding when and where to point a sensor to re-acquire a tracked object\n",
"\n",
"Under the hood, satkit propagates the 6x6 state transition matrix $\\Phi(t, t_0)$ alongside the equations of motion. The covariance at time $t$ is then computed as:\n",
"\n",
"$$P(t) = \\Phi(t, t_0) \\, P(t_0) \\, \\Phi(t, t_0)^T$$\n",
"\n",
"where $P(t_0)$ is the initial covariance matrix."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-04-05T17:29:36.373076Z",
"iopub.status.busy": "2026-04-05T17:29:36.372958Z",
"iopub.status.idle": "2026-04-05T17:29:36.595333Z",
"shell.execute_reply": "2026-04-05T17:29:36.595112Z"
}
},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import math\n",
"import matplotlib.pyplot as plt\n",
"import scienceplots # noqa: F401\n",
"plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
"%config InlineBackend.figure_formats = ['svg']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setting Up a Satellite State\n",
"\n",
"We create a satellite in a near-circular LEO orbit at roughly 400 km altitude. The position and velocity are specified in the GCRF (Geocentric Celestial Reference Frame) inertial frame."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-04-05T17:29:36.596599Z",
"iopub.status.busy": "2026-04-05T17:29:36.596493Z",
"iopub.status.idle": "2026-04-05T17:29:36.598931Z",
"shell.execute_reply": "2026-04-05T17:29:36.598709Z"
}
},
"outputs": [],
"source": [
"# Epoch\n",
"t0 = sk.time(2024, 6, 15, 12, 0, 0)\n",
"\n",
"# LEO orbit: ~400 km altitude, near-circular\n",
"r = 6.778e6 # meters (Earth radius + 400 km)\n",
"v = math.sqrt(sk.consts.mu_earth / r) # circular orbit speed\n",
"\n",
"pos = np.array([r, 0.0, 0.0]) # meters, GCRF\n",
"vel = np.array([0.0, v, 0.0]) # m/s, GCRF\n",
"\n",
"state = sk.satstate(t0, pos, vel)\n",
"print(f\"Position: {state.pos} m\")\n",
"print(f\"Velocity: {state.vel} m/s\")\n",
"print(f\"Altitude: {np.linalg.norm(state.pos) - 6.378e6:.0f} m\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setting Initial Covariance\n",
"\n",
"Uncertainty is most naturally expressed in a satellite-local frame where the axes are aligned with the orbit geometry. `satkit` provides a unified API for this — `set_pos_uncertainty(sigma, frame)` and `set_vel_uncertainty(sigma, frame)` — that accepts any of the supported orbital or inertial frames: `LVLH`, `RIC` (a.k.a. RSW/RTN), `NTW`, or `GCRF`.\n",
"\n",
"Here we use **LVLH** (Local Vertical, Local Horizontal), where z = nadir, y = anti-angular-momentum, and x ≈ along-track for a circular orbit. Calling `set_pos_uncertainty` and `set_vel_uncertainty` in sequence builds up the full 6×6 covariance — each call preserves the block it isn't updating."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-04-05T17:29:36.614467Z",
"iopub.status.busy": "2026-04-05T17:29:36.614355Z",
"iopub.status.idle": "2026-04-05T17:29:36.616555Z",
"shell.execute_reply": "2026-04-05T17:29:36.616348Z"
}
},
"outputs": [],
"source": [
"# 1-sigma position uncertainty in LVLH (x ≈ along-track, y = cross-track, z = nadir)\n",
"sigma_pos_lvlh = np.array([50.0, 20.0, 10.0])\n",
"state.set_pos_uncertainty(sigma_pos_lvlh, frame=sk.frame.LVLH)\n",
"\n",
"# 1-sigma velocity uncertainty in LVLH\n",
"sigma_vel_lvlh = np.array([0.05, 0.02, 0.01]) # m/s\n",
"state.set_vel_uncertainty(sigma_vel_lvlh, frame=sk.frame.LVLH)\n",
"\n",
"print(\"Position covariance block (GCRF, stored internally), meters^2:\")\n",
"print(state.cov[0:3, 0:3])\n",
"print(\"\\nFull 6x6 covariance matrix (GCRF):\")\n",
"print(state.cov)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Propagating with Covariance\n",
"\n",
"When a `satstate` has a covariance matrix set, calling `propagate` will propagate both the orbit and the covariance using the state transition matrix. Let's propagate forward by one orbital period and examine how the covariance changes."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-04-05T17:29:36.617523Z",
"iopub.status.busy": "2026-04-05T17:29:36.617442Z",
"iopub.status.idle": "2026-04-05T17:29:36.700065Z",
"shell.execute_reply": "2026-04-05T17:29:36.699833Z"
}
},
"outputs": [],
"source": [
"# Orbital period\n",
"a = np.linalg.norm(state.pos) # semi-major axis (circular orbit)\n",
"period_sec = 2.0 * math.pi * math.sqrt(a**3 / sk.consts.mu_earth)\n",
"print(f\"Orbital period: {period_sec / 60:.1f} minutes\")\n",
"\n",
"# Propagate forward by one orbit\n",
"state_1orbit = state.propagate(sk.duration.from_seconds(period_sec))\n",
"\n",
"print(\"\\nInitial 1-sigma position uncertainty (GCRF diagonal), meters:\")\n",
"print(f\" {np.sqrt(np.diag(state.cov[0:3, 0:3]))}\")\n",
"\n",
"print(\"\\nAfter 1 orbit, 1-sigma position uncertainty (GCRF diagonal), meters:\")\n",
"print(f\" {np.sqrt(np.diag(state_1orbit.cov[0:3, 0:3]))}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Examining Uncertainty Growth Over Time\n",
"\n",
"To understand how uncertainty evolves, we propagate the state at regular intervals over 24 hours and extract the total position uncertainty (RSS of the 1-sigma values) at each step."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-04-05T17:29:36.701069Z",
"iopub.status.busy": "2026-04-05T17:29:36.700997Z",
"iopub.status.idle": "2026-04-05T17:29:37.044898Z",
"shell.execute_reply": "2026-04-05T17:29:37.044601Z"
}
},
"outputs": [],
"source": [
"# Propagate at 5-minute intervals over 24 hours\n",
"dt_minutes = 5\n",
"num_steps = int(24 * 60 / dt_minutes)\n",
"\n",
"times_hr = np.zeros(num_steps)\n",
"sigma_pos_rss = np.zeros(num_steps) # RSS position uncertainty\n",
"sigma_x = np.zeros(num_steps)\n",
"sigma_y = np.zeros(num_steps)\n",
"sigma_z = np.zeros(num_steps)\n",
"\n",
"for i in range(num_steps):\n",
" dt = sk.duration.from_minutes(dt_minutes * (i + 1))\n",
" s = state.propagate(dt)\n",
" \n",
" times_hr[i] = dt_minutes * (i + 1) / 60.0\n",
" pos_cov = s.cov[0:3, 0:3]\n",
" sigmas = np.sqrt(np.diag(pos_cov))\n",
" sigma_x[i] = sigmas[0]\n",
" sigma_y[i] = sigmas[1]\n",
" sigma_z[i] = sigmas[2]\n",
" sigma_pos_rss[i] = np.sqrt(np.trace(pos_cov))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-04-05T17:29:37.046173Z",
"iopub.status.busy": "2026-04-05T17:29:37.046096Z",
"iopub.status.idle": "2026-04-05T17:29:37.238067Z",
"shell.execute_reply": "2026-04-05T17:29:37.237834Z"
}
},
"outputs": [],
"source": [
"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)\n",
"\n",
"# Total RSS position uncertainty\n",
"ax1.plot(times_hr, sigma_pos_rss / 1e3, linewidth=1.5)\n",
"ax1.set_ylabel(\"RSS Position 1-sigma (km)\")\n",
"ax1.set_title(\"Position Uncertainty Growth Over 24 Hours\")\n",
"ax1.grid(True, alpha=0.3)\n",
"\n",
"# Per-component uncertainty (GCRF)\n",
"ax2.plot(times_hr, sigma_x / 1e3, label=\"X (GCRF)\", linewidth=1.0)\n",
"ax2.plot(times_hr, sigma_y / 1e3, label=\"Y (GCRF)\", linewidth=1.0)\n",
"ax2.plot(times_hr, sigma_z / 1e3, label=\"Z (GCRF)\", linewidth=1.0)\n",
"ax2.set_xlabel(\"Time (hours)\")\n",
"ax2.set_ylabel(\"Position 1-sigma (km)\")\n",
"ax2.set_title(\"Per-Component Position Uncertainty (GCRF)\")\n",
"ax2.legend()\n",
"ax2.grid(True, alpha=0.3)\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The LVLH Frame\n",
"\n",
"The LVLH (Local Vertical, Local Horizontal) frame is a coordinate system attached to the orbiting satellite, defined as:\n",
"\n",
"- **x-axis**: approximately along-track (in the direction of velocity)\n",
"- **y-axis**: cross-track (opposite to the orbital angular momentum vector $\\mathbf{h} = \\mathbf{r} \\times \\mathbf{v}$)\n",
"- **z-axis**: nadir (pointing toward Earth center, i.e., $-\\hat{r}$)\n",
"\n",
"This frame is useful because orbital uncertainties have a natural physical interpretation in LVLH: along-track errors tend to grow fastest (due to period uncertainty), while radial errors remain relatively bounded.\n",
"\n",
"The `qgcrf2lvlh` property on `satstate` provides the quaternion that rotates vectors from the GCRF inertial frame into the LVLH frame."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-04-05T17:29:37.239161Z",
"iopub.status.busy": "2026-04-05T17:29:37.239096Z",
"iopub.status.idle": "2026-04-05T17:29:37.241130Z",
"shell.execute_reply": "2026-04-05T17:29:37.240910Z"
}
},
"outputs": [],
"source": [
"# Demonstrate the LVLH quaternion\n",
"q_lvlh = state.qgcrf2lvlh\n",
"print(\"Quaternion GCRF -> LVLH:\", q_lvlh)\n",
"\n",
"# Verify the LVLH axes:\n",
"# z-axis should point toward nadir (-r_hat)\n",
"r_hat = state.pos / np.linalg.norm(state.pos)\n",
"z_lvlh = q_lvlh * r_hat # Should be [0, 0, -1]\n",
"print(f\"\\nr_hat in LVLH (expect ~[0, 0, -1]): {z_lvlh}\")\n",
"\n",
"# y-axis should be along -h_hat\n",
"h = np.cross(state.pos, state.vel)\n",
"h_hat = h / np.linalg.norm(h)\n",
"y_lvlh = q_lvlh * h_hat # Should be [0, -1, 0]\n",
"print(f\"h_hat in LVLH (expect ~[0, -1, 0]): {y_lvlh}\")\n",
"\n",
"# x-axis should be approximately along velocity\n",
"v_hat = state.vel / np.linalg.norm(state.vel)\n",
"x_lvlh = q_lvlh * v_hat # Should be ~[1, 0, 0]\n",
"print(f\"v_hat in LVLH (expect ~[1, 0, 0]): {x_lvlh}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Covariance in the LVLH Frame\n",
"\n",
"While satkit stores and propagates the covariance in GCRF, we can rotate it into the LVLH frame at any point to examine the uncertainty in physically meaningful orbit-relative directions."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-04-05T17:29:37.242073Z",
"iopub.status.busy": "2026-04-05T17:29:37.242015Z",
"iopub.status.idle": "2026-04-05T17:29:37.592986Z",
"shell.execute_reply": "2026-04-05T17:29:37.592658Z"
}
},
"outputs": [],
"source": [
"def cov_gcrf_to_lvlh(sat):\n",
" \"\"\"Rotate the 3x3 position covariance from GCRF to LVLH.\"\"\"\n",
" q = sat.qgcrf2lvlh\n",
" # Build rotation matrix by rotating basis vectors\n",
" R = np.column_stack([\n",
" q * np.array([1.0, 0.0, 0.0]),\n",
" q * np.array([0.0, 1.0, 0.0]),\n",
" q * np.array([0.0, 0.0, 1.0]),\n",
" ])\n",
" pos_cov_gcrf = sat.cov[0:3, 0:3]\n",
" return R @ pos_cov_gcrf @ R.T\n",
"\n",
"\n",
"# Propagate and extract LVLH uncertainties\n",
"sigma_along = np.zeros(num_steps)\n",
"sigma_cross = np.zeros(num_steps)\n",
"sigma_radial = np.zeros(num_steps)\n",
"\n",
"for i in range(num_steps):\n",
" dt = sk.duration.from_minutes(dt_minutes * (i + 1))\n",
" s = state.propagate(dt)\n",
" cov_lvlh = cov_gcrf_to_lvlh(s)\n",
" sigmas = np.sqrt(np.diag(cov_lvlh))\n",
" sigma_along[i] = sigmas[0] # x = along-track\n",
" sigma_cross[i] = sigmas[1] # y = cross-track\n",
" sigma_radial[i] = sigmas[2] # z = radial"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-04-05T17:29:37.594181Z",
"iopub.status.busy": "2026-04-05T17:29:37.594104Z",
"iopub.status.idle": "2026-04-05T17:29:37.692983Z",
"shell.execute_reply": "2026-04-05T17:29:37.692760Z"
}
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(10, 5))\n",
"\n",
"ax.plot(times_hr, sigma_along / 1e3, label=\"Along-track\", linewidth=1.5)\n",
"ax.plot(times_hr, sigma_cross / 1e3, label=\"Cross-track\", linewidth=1.5)\n",
"ax.plot(times_hr, sigma_radial / 1e3, label=\"Radial\", linewidth=1.5)\n",
"\n",
"ax.set_xlabel(\"Time (hours)\")\n",
"ax.set_ylabel(\"Position 1-sigma (km)\")\n",
"ax.set_title(\"Position Uncertainty in LVLH Frame Over 24 Hours\")\n",
"ax.legend()\n",
"ax.grid(True, alpha=0.3)\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The LVLH-frame plot reveals the characteristic pattern of orbital uncertainty growth: **along-track uncertainty dominates** and grows roughly linearly with time, because a small error in velocity (and hence orbital period) causes the satellite to drift ahead or behind its predicted position. Cross-track and radial uncertainties remain comparatively bounded."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}